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A new spectral type method for solving the one dimensional quantum- 
mechanical Lippmann-Schwinger integral equation in configuration space is 
described. The radial interval is divided into partitions, not necessarily of 
equal length. Two independent local solutions of the integral equation are 
obtained in each interval via Clenshaw-Curtis quadrature in terms of Cheby- 
shev Polynomials. The local solutions are then combined into a global solution 
by solving a matrix equation for the coefficients. This matrix is sparse and the 
equation is easily soluble. The method shows excellent numerical stability, as 
is demonstrated by several numerical examples. 

I. INTRODUCTION. 

In conjunction with Professor I. Koltracht and several students in the Mathematics De- 
partment of the University of Connecticut, we are in the process of developing a method 

for solving the one- dimensional Lippmann-Schwinger integral equations in configuration 
space, associated with the corresponding Schrodinger equation. At a later stage we hope 
to generalize the method to other types of integral equations, and attempt to increase the 
dimension of the number of variables. Our method is excellently suited for cases where either 
the potential has a long range with many oscillations occurring in the wave function, or for 
large systems of coupled equations in which some of the channel energies are positive and 
others negative. In both situations the conventional finite difference methods (Numerov, for 
example) experience difficulties because of the larger accumulation of roundoff errors, and 
because of linear independence problems in implementing the asymptotic boundary condi- 
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tions in the coupled channel case. A good physical application of our method is very likely 
the process of photoassociation in the collision of cold atoms, because of the very long range 
of the interaction between the atoms, and also because of the many excited configurations 
which participate and which lead to many coupled channels. Photoassociation is a relatively 
new field of study. It was recently reviewed by Julienne 0, and is observed experimen- 
tally at various laboratories, including the University of Connecticut ||. The calculation of 
three-body reactions in coordinate space might furnish another example because of the large 
distances involved, especially if we are able to generalize our method to two dimensions, but 
this conjecture has not yet been examined. 

The Chebyshev expansion method is "spectral", i.e., the accuracy increases faster than 
any inverse power of the number of mesh points employed. Our method is based on that 
of Greengard and Rokhlin but differs |IJ from it substantially in the manner in which 
the local solutions in each partition are combined into the global solution. We prefer to 
use Chebyshev polynomials for the expansion basis in each partition because these functions 
have excellent integral properties || , and the node points can be obtained in terms of simple 
algebraic cosine expressions, rather than in terms of non-analytic expressions required for 
the zeros of Legendre Polynomials, for example. 

II. THE ALGORITHM. 

For the one-channel case our integral equation method (IEM) is as follows. The 
Schrodinger equation for the partial wave function ip{R) 

+ V L (R) - fc 2 ) MR) = (1) 

is first transformed into an equivalent integral equation 

iJ>(R) = F(R) + [ g (R, R')V L {R')i>(R') dR'. (2) 
Jo 

In configuration space the Green's function Go is semi-separable, i.e., it is given by the 
product of two independent solutions F and G of the Schrodinger equation, i.e., Go(R, R') oc 
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Fo(R<) x Go(R < ). Here R < (R>) is the smaller (larger) of R and R', F = sin(Ar), Go = 
cos(kr), and F(R) is the driving function which is equal to F in the uncoupled case, and 
the upper limit in the integral is T = Rmox- Physicists would be tempted to use instead of 
F Q and Go the L-dependent Riccati-Bessel functions Fl and Gl, but because the spectral 
expansion technique is very robust, it allows us to place the L(L + l)/R 2 term into the 
potential Vl 

V L (R)=L(L + 1)/R 2 + V(R), (3) 

without loss of accuracy [fl]. 

We avoid the occurrence of large non-sparse matrices by: a) Dividing the integration 
interval [0, Rmclx] into m partitions [0, bi], b 2 ], ... [6 m -i, Rmox]- The size of each partition 
can be arbitrary, but two or three partitions per local wave number is optimum, b) Solving 
the integral equation separately in each partition i for two functions Y and Z 

Yi- 1* So{R, R')V L {R')Y t {R')dR' = F(R), <R<h (4) 

Zi- r G (R, R')Vl(R') Zi (R') dR' = G(R), fc-i < R < h. (5) 

Jbi-i 

The method of solution for this step leads to small non-sparse matrices. It consists in 
expanding the unknown solutions into Chebyshev polynomials and making use of their 
convenient properties [4-6]. c) "Stitching" together the separate solutions into the global 
one by means of the expressions, valid in each partition 

ip(R) = AiYi(R) + BiZi(R), h^KRKh. (6) 

This " stitching" procedure leads to a large but sparse matrix. It produces a seamless smooth 
continuation of the function if} from one partition into the next, and is a consequence of the 
semi-separability of Q. As a consequence of this property the coefficients A and B obey 
the matrix equation 
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where the ctj and w are two-rowed columns 
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For the one-channel case the Mjj are 2x2 matrices whose elements are formed out 
of overlap integrals between combinations of the functions Y, Z, the potential V and the 
functions F and G We call the matrix formed out of the blocks of the M^s the "Big 
Matrix". It is sparse, and of dimension 2m x 2m, where m is the number of partitions, 
while the matrices required to solve the equations for the Y and Z in each partition are 
of the dimension n x n, where n is the number of Chebyshev points in each partition. 
Following standard practice, that number is chosen to have the value n = 16. The end 
result of the whole calculation is the set of coefficients of the Chebyshev polynomials for the 
functions Y and Z in each partition, as well as the coefficients A and B for each partition. 
The calculations in one partition of the Y, Z, and the overlap integrals required for the 
construction of the My can be performed independently of each other. Hence this part of 
the calculation can be carried out in a parallel computer architecture. Because of the simple 
structure of the "Big Matrix", the solution of Eq. (7) can be performed by simple pivoting, 
and does not require special big matrix techniques. The values of ip and its derivative 
dijj/dR can be obtained at any arbitrary point from the solution given by Eq. (6), because 
the values and the derivatives of Chebyshev functions are known analytically. At the end 
point Ruax — T the expression for dip/dR is especially simple. 

When the Schrodinger Eq. is augmented to a system of coupled equations with N 



channels, then each Fq and the corresponding Y is augmented into N column vectors of 
length N, and similarly for the G's and Z's. The size of the block matrices My increases 
correspondingly to the dimension 2N x 2N, but the sparse structure of the "Big Matrix" 
remains the same. For the channels which are closed, the functions Fq and Go are replaced by 
sinh(Ki?) and exp(— kR), respectively, where k is the imaginary part of the asymptotic wave 
number in the closed channel. Numerical experiments 0] for L = show that the stability 
and high precision of the results is maintained in this case also, provided that a scaling 
procedure is introduced which reduces the numerical disparity between the two exponentials 
exp(±/?i?) for large values of kR. If in the positive energy channels the angular momentum 
L 7^ 0, it is still advantageous to use for the channel Green functions the undistorted Q ones, 
even though one now has to solve the coupled integral equations as many times as there are 
open channels, each time with a different driving function F. The big matrix M is the same 
in each case, the only change is in the vector u on the right hand side of Eq. (7), and one 
can show |7j that the corresponding solutions are linearly independent asymptotically. In 
the next section several examples illustrate the numerical properties of our method. 



III. ACCUMULATION OF ROUNDOFF ERRORS. 

This property is demonstrated by a numerical example for an uncoupled channel case 
with L ^ 0, but in which the only potential present is L(L + 1)/R 2 . In this case the solution 
is a Riccati-Bessel function Fl(R) = kRiiikR). The numerical solution ip of Eq. (2) is 
proportional to Fl, and the constant of proportionality is obtained from the Wronskian of 
ip with Gi at R — Rmcix- The calculations are done in double precision on a IBM 3090 
Mainframe machine. The functions Gl and Fl at Ruax are obtained from the International 
Mathematical Scientific Library (IMSL). The error in the numerical solution ip is obtained 
by comparing it with the IMSL values of F^ at many points in the interval [0, Rmo,x], and 
the maximum of the absolute values of all these differences is denoted as "Error". This 
error is plotted in Fig. 1 as a function of N, the total number of integration points in 
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[0, Rmclx}- For the IEM that number is equal to 16 times the number of partitions m, while 
for the Numerov method it is Ruax divided by the uniform mesh size. In this example 
with k = l/m -1 , Ruax = 50 fm and L = 6, one sees that the error rapidly drops as a 
function of N, and reaches a minimum (equal to machine accuracy in this case) at the point 
where the roundoff error overtakes the truncation error. Beyond that minimum the rise is a 
measure of the rapidity of accumulation of roundoff error. One sees from the figure that the 
Numerov error decreases much more slowly with N, and the accumulation of round-off error 
increases much faster. This is the reason why the best accuracy for the Numerov method 
is much worse. A variable step size improved Numerov method due to Raptis and Cash || 
accumulates roundoff error more slowly than the Numerov method, but still does not reach 
the quality of the IEM. 

A more taxing example with L = 8, k = 40fm~ 1 and Rmox — 50 (not illustrated here) 
shows a similarly slow accumulation of the roundoff error for the IEM. In this case there 
are 580 nodes in the wave function, and the smallest error is one order of magnitude larger 
than machine accuracy JTJ. The corresponding value of N is 12800, which corresponds to 
800 partitions in the [0, Rmux] interval, or 20 Chebyshev points per half wavelength. The 
corresponding number for best accuracy (~ 10~ 8 ) of the Numerov method is 640 points per 
half wavelength, larger by a factor of 30. There results are summarized in the Table below. 



Maximum Accuracy of the Riccati-Bessel Function. 







# oscill. 


# of points/(local A) 


k(fm- x ) 


KM 


0<R<50 


IEM Numerov 


1 


6.3 


8 


50 1260 


40 


0.16 


310 


40 1280 



IV. IMPLEMENTATION OF BOUNDARY CONDITIONS. 

For the case of one channel the asymptotic boundary conditions ip(R) ~ F l (R)+wGl(R) 
are easily obtained in both the IEM and the finite difference methods (FDM), by simple 
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normalization at the end point. However if there are two or more coupled channels, then 
the various supposedly independent wave functions obtained in the FDM by starting the 
solutions near the origin independently, can loose a large part of their independence near 
the matching point, and the process of obtaining the appropriate linear combination which 
satisfies the desired boundary condition looses accuracy. This is not the case for the IEM, 
as will now be demonstrated for the case L = 0. 

In this example @ we have two channels, both with the same positive energy 
E = k 2 fm~ 2 , and all potentials (coupling as well as diagonal) are of exponential form 
Voexp(— r/a), with the same decay parameter a and the same magnitude of the strength Vq. 
For the single channel L = case there exists an analytic solution given in terms of Bessel 
functions J u (y) of imaginary index v = ±2iak, and argument y = 2ay/—Vo exp(— r/2a), 
which can be generalized to the coupled channel case when the energies in all the channels 
are equal. In our example | Vq |= 5/\/2fm~ 2 and a = 4 fin. Channel 1 is the incident 
channel, with V > 0, and ipi(R) ~ F (R) + wiGq(R), while in channel 2 V < and 
ip2{R) ~ WzGo(R). The value of Ruax is 180/m and 140/m, respectively, in the IEM and 
Numerov calculations, again carried out in double precision on the IBM mainframe. For 
each value of the wave number k the number of points in each method of calculation is 
varied until the combined error in w\ and u>2 is a minimum, and this minimum error is 
displayed on the vertical axis in Fig. 2. The reason why the Numerov error (open circles) 
increases as k decreases is because of the lack of independence of the two solutions. This is 
shown either by an analytical argument f?|, valid in this exponential case, or can seen from 
the fact that the Numerov error for w in the uncoupled channels is much less dependent on 
k. The latter is shown by the lines with squares or triangles, labeled "Attr." and "Rep.", 
respectively. The error in the IEM is even more independent on k, and is also much smaller 
since the boundary conditions are automatically built into the Green's functions. 

The case with the same exponential potentials in which however the energy in the second 
channel is negative has also been examined. Since an analytical solution does not exist in 
this case, the accuracy of the solution has to be inferred from the stability of the values of 
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w\ and W2- For the IEM the values of both w\ and W2 are stable to 13 sign figures, while 
for the Numerov method the value of W2 is correct to less than 4 significant figures for 
kRuax > 25, .while W\ is correct to about 7. 

V. SUMMARY AND CONCLUSIONS. 

We have demonstrated a method for solving the Lippmann-Schwinger integral equation 
in configuration space which is linear in the number of the integration points, and which 
is numerically very stable. The stability of the IEM method is due to two factors: 1. The 
solution of integral equations lead to a smaller accumulation of roundoff errors than the 
finite difference methods of solving a differential equation, and 2), The IEM is a spectral 
method, which has an inherently higher accuracy than finite difference methods. 

Furthermore the IEM is easy to implement: a) The boundary conditions on the channel 
wave functions are automatically incorporated through the Green's functions; b) The size 
of each partition — fej can be assigned independently of the sizes of the other partitions; 
c) The overlap integrals required for the construction of matrix elements can be obtained 
with the Gauss-Legendre quadrature because the functions in the integrand, being given 
in terms of Chebyshev expansions, can be calculated easily for any given sets of points, 
no interpolation being required. The Curtis-Clenshaw (C-C) quadrature || is also a good 
option (see Table 2 in Ref. because the wave functions and potentials are already known 
at the appropriate Chebyshev points, and certainly exceeds the efficiency of equidistant 
point methods. 
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B. Figures 

Fig. 1 Accumulation of roundoff 
errors for a Riccati-Bessel 
function, as a function of 
the number of integration steps N. 

Fig. 2. Accuracy of the R-matrix 
elements for two channels 
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coupled by exponential potentials, 
as a function of the wave number. 



